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THE  IMPORTANCE  OF  BOUNDARY  CONDITIONS 


IN  THE  NUMERICAL  TREATMENT  OF  HYPERBOLIC  EQUATIONS 

Gino  Moretti 

Polytechnic  Institute  of  Brooklyn 
Farmingdale,  New  York 

ABSTRACT 

Many  of  the  existing  computations  of  initial-and- 
boundary  value  problems  in  fluid  mechanics  suffer  from 
unrealistic  treatment  of  boundary  points.  Three  categories 
of  boundaries  are  briefly  discussed:  rigid  walls,  arbitrary 
boundaries  of  a  computational  region  in  a  subsonic  flow, 
and  shock  waves.  An  attempt  is  made  to  show  in  what 
sense  the  numerical  treatment  of  such  boundaries  may  be 
physically  wrcng  and  what  can  be  done  instead.  Examples 
from  the  blunt  body  problem,  the  transonic  flow  in  a  nozzle, 
the  incompressible  inviscid  flow  past  a  circle,  and  the  quasi- 
one-dimensional  flow  in  a  Laval  nozzle,  are  shown. 

r . 


This  paper  was  presented  at  the  International  Symposium 
on  High-Speed  Computing  in  Fluid  Dynamics,  Monterey, 
California,  August  18-24,  1968. 


1.  INTRODUCTION 


A  great  deal  of  work  has  been  performed  in  the 
attempt  to  use  numerical  techniques  to  solve  complicated 
problems  in  fluid  mechanics,  but  too  many  results  are  far 
from  being  satisfactory.  Failures  are  sometimes  ascribed 
to  such  mysterious  causes  as  "non-linear  instability"  (a 
term  which  is  meaningless  for  lack  of  a  definition  and  sounds 
very  much  like  the  "Hie  sunt  leones"  label  attached  to  un¬ 
explored  lands  in  middle-age  maps).  I  will  try  to  show  in 
this  paper  that  in  mixed  initial -and- boundary -value  prob¬ 
lems  major  troubles  arise  if  the  boundary  conditions  are  not 
properly  handled. 

Surprisingly,  the  difficulties  of  boundary -value  prob¬ 
lems  so  far  seem  to  have  been  overestimated.  Let  us  quote, 
for  example,  from  page  128  of  the  authoritative  book  by  Richt- 
myer  and  Morton*:  "For  one  thing  the  effect  of  boundary  con¬ 
ditions  has  not  so  far  been  considered";  and  from  page  167; 
"To  treat  the  more  complex  boundaries  that  can  arise  with 
more  space  dimensions  is  very  much  more  difficult  and  little 
has  so  far  been  achieved  in  this  area".  Consequently,  Richt- 
myer's  book,  which  is  a  study  in  mathematical  foundations, 
does  not  even  attempt  to  attack  the  problem.  Unfortunately, 
since  nobody  else  does  it  in  a  book  of  the  same  level,  oriented 
towards  applications,  the  physicist  and  the  engineer  who  face 
the  task  of  solving  a  practical  problem  are  led  to  believe  that 
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what  is  not  discussed  by  Richtmyer  is  a  matter  of  no  con¬ 
sequence.  In  other  words,  they  tend  to  unde  res  timate  the 
importance  of  boundary  conditions  which  are,  instead,  the 
governing  elements  of  the  entire  computation. 

The  situation  is  worsened  by  the  fact  that  the  few 
examples  of  multi-dimensional  problems  given  by  Richt¬ 
myer  in  Chapter  13  (taken  from  the  existing  literature) 
show  a  number  of  troubles  within  the  computational  region 
whereas  such  troubles  are  generated  at  the  boundaries. 

There  are,  in  fact,  two  classes  of  points  to  con¬ 
sider  in  a  numerical  computation;  interior  points  and 
boundary  points.  As  far  as  the  interior  points  are  con¬ 
cerned,  many  techniques  are  available.  One  can  discuss 
their  relative  advantages  or  disadvantages  in  terms  of,  say, 
computational  speed,  accuracy,  programming  difficulty 
and  iio  on.  There  are  ways  of  determining  such  properties 
through  a  local  analysis  and  some  of  the  techniques,  accurate 
to  the  second  order,  prove  to  be  very  good  for  all  practical 
purposes.  Nevertheless,  many  of  the  published  computations 
show  a  progressive  deterioration  of  results  which  is  incon¬ 
sistent  with  the  good  local  qualities  of  the  technique  used  for 
interior  points  and  is  then  attributed  to  "non-linear  insta¬ 
bility".  A  "cure"  for  the  latte*  is  sought  in  damping  devices 
("artificial  viscosity"),'  a  procedure  which  is  physically  un¬ 
justified  and  which  may  provide  smooth  results  at  a  very  high 
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price,  that  is,  the  destruction  of  accuracy  .  In  what  follows, 
I  would  like  to  show  a  number  of  cases  where  there  is  no  rea¬ 
son  for  introducing  a  concept  of  non-linear  instability,  con¬ 
flicting  with  the  local  stability,  if  care  is  taken  not  to  feed 
unrealistic  perturbations  from  the  boundaries.  To  empha¬ 
size  the  relevant  arguments,  the  presentation  will  necessarily 
be  sketchy  and  many  details  will  be  left  unmentioned. 

2.  BOUNDARY  CONDITIONS  ON  A  RIGID  WALL 

Let  us  begin  witn  few  obvious  statements. 

(1)  A  fluid  flow  problem  is  described  by  a  system  of 
partial  differential  equations  which  are  classically  known  as 
the  indefinite  equations  of  motion.  Such  equations  are  called 
indefinite  because  they  apply  to  any  problem  in  general  but 
the  equations  per  se  do  not  define  a  specific  problem. 

(2)  The  problem  is  defined  only  when  a  proper  set  of 
initial  and/or  boundary  conditions  is  given.  The  boundary 
conditions  are  such  an  important  part  of  a  definition  of  a 
problem  that  the  patterns  of  two  flow  fields  can  be  completely 
different  from  one  another  because  of  some  differences  in 
the  flow  boundaries,  despite  the  fact  that  both  flows  obey  the 
same  system  of  indefinite  partial  differential  equations. 

(3)  For  each  system  of  equations  there  are  a  number 
of  necessary  and  sufficient  boundary  conditions.  For  example, 
for  an  inviscid  flow  the  condition  on  a  rigid,  fixed  wall  is  the 
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vanishing  of  the  normal  component  of  the  velocity,  V  . 

At  the  risk  of  appearing  trivial,  let  me  state  that  no 
other  condition  may  be  Imposed  on  the  rigid  wall  since  the  one 
above  is  sufficient  to  determine  the  flow  field. 

Now,  let  us  see  what  happens  when  a  problem  of  in- 
viscid  flow  is  treated  by  a  numerical  technique,  A  certain 
mesh  is  used.  At  the  interior  mesh-points  the  partial  differ¬ 
ential  equations  are  substituted  by  finite -difference  equations 
(and  we  will  assume  that  the  substitution  satisfies  all  conditions 
for  local  stability  and  convergence).  At  each  computational 
step,  information  is  transmitted  from  each  point  to  its  neigh¬ 
boring  points  via  the  computation  of  finite  differences.  In  this 
way,  the  boundary  mesh-points  influence  their  neighbor?  and 
transmit  the  effects  of  the  boundary  condition  into  the  flow 
field. 

For  this  to  happen,  the  values  of  all  physical  quantities 
at  boundary  points  must  be  known  at  each  computational  step. 

At  a  rigid,  fixed  boundary  only  the  normal  component  of  the 
velocity,  is  clearly  defined,  but  the  tangential  velocity, 

V^,  the  pressure,  p,  the  density,  P,  are  also  required.  At 
this  point,  we  are  required  to  find  some  way  of  computing 
them  using  information  from  interior  points  plus  the  boundary 
condition.  In  other  words,  we  must  solve  finite  difference 
approximations  to  the  partial  differential  equations  of  motion 
in  their  limiting  form  as  one  approaches  the  boundary.  If 
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we  do  not  proceed  properly,  we  face  the  risk  of  over- specify¬ 
ing  the  boundary  conditions  themselves  and,  in  all  probability, 
these  overspecified  boundary  conditions  will  not  be  consistent 
with  the  nature  of  the  boundary  and  the  limiting  forms  of  the 
equations  of  motion. 

Unexpected  procedures  are  found  in  the  literature. 

In  general,  it  seems  that  an  "easy"  way  for  solving  the  equa¬ 
tions  of  motion  at  points  on  a  rigid  wall  consists  of  adding  an 
extra  row  of  points  behind  the  wall  and  computing  the  wall 
points  with  the  same  routine  used  for  interior  points.  At 
each  computational  step,  the  problem  is  then  shifted  into 
that  of  defining  values  at  the  extra  points  which  provide  good 
derivatives  at  the  wall  points.  To  this  effect,  many  authors 
use  what  they  call  a  "reflection-principle"  which  iff  not  a 
principle  at  all  and  should  rather  be  called  a  reflection  tech¬ 
nique.  I  have  been  unable  to  find  a  precise  definition  of  such 
a  technique  in  the  literature.  In  its  crudest  formulation,  all 
physical  parameters  are  specularly  reflected  on  the  rigid  wall 

except  which  is  assumed  as  anti -symmetric.  Bohachevsky 

3  4 

et  al.  ,  page  603  and  605;  Bohachevsky  et  al.  ,  page  778; 

5  6 

Bur  stein  ,  page  2114;  Singleton  ,  page  2-2).  Unfortunately, 
only  the  assumption  regarding  Vn  is  legitimate  since  it 
produces  a  vanishing  Vn  at  the  wall  automatically  and  the 
normal  derivatives  of  Vn  so  computed  is  the  same  as  a 
one-sided  derivative  computed  by  finite  differences,  using 
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the  interior  point  information  (the  accuracy  so  obtained  is, 
however,  very  low).  The  other  assumptions  force  the 
normal  derivatives  of  the  remaining  quantities  to  vanish 
at  the  wall.  These  are  redundant  conditions  and  are 
physically  wrong. 

To  show  how  far  the  real  situation  is  from  the 
assumed  one,  let  us  consider  the  case  of  a  circular  ob¬ 
stacle  in  a  uniform  flow,  the  flow  being  steady,  two- 
dimensional,  incompressible,  and  irrotational.  With  a 
proper  scaling  of  the  variables,  the  flow  is  described  by 
the  complex  potential, 

W  =  z  +  I 
z 

where  z  =  x  +  iy  is  the  complex  coordinate  in  the  physical 
plane.  Consequently,  the  complex  velocity  is 


The  velocity  components  in  polar  coordinates  (u  in  the 
radial  direction,  v  in  the  transverse  direction)  are 

US  (1  ■  i.  )  COS  0 
r2 

v  =  (  1  +  I  )  sin  0 
r2 

and  the  pressure  coefficient  is 

C  1  +  2  cos  2  9 
P  "  "r“  ra 
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The  obstacle  is  defined  by  r  =  1.  The  normal  derivatives 
of  pressure  and  velocity  component  at  tin  Jtacle  wall 
are: 


=  8  sin3  0 
Sr 


Su 

Sr 


=  2  cos  6 


Sv 

Sr 


=  -  2  sin  0 


at  r  =  1 


Obviously,  such  normal  derivatives  of  pressure  and  v  are 
far  from  vanishing,  except  at  the  stagnation  points  (0  =  0,  tt)  . 
For  example,  Fig.  1  shows  the  distribution  of  Cp  and  v  along 
the  radius  at  0  =  tt/2  (solid  lines). 

The  flow  field  described  above  can  be  considered  as  the 
asymptotic  state  of  a  time -dependent  computation.  If,  in 
performing  the  latter,  conditions  are  imposed  at  the  bound¬ 
ary  which  are  patently  conflicting  with  the  real  nature  of 
the  solution,  one  is  left  to  wonder  why  certain  computations 
did  not  "explode".  There  are  two  reasons.  One  is  that  the 
departure  from  symmetry  is  not  substantial  in  the  vicinity 
of  a  stagnation  point  or  along  a  wall  which  does  not  differ 
too  much  from  an  infinite  straight  line.  The  other,  and 
more  important,  is  that  a  very  high  artificial  viscosity  has 
been  used  --  an  ill-advised  procedure  as  I  said  before. 
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Results  of  the  following  example  support  the  above 

arguments.  Consider  a  blunt  body  with  a  circular  nose  in 

a  supersonic  stream  (M^  =  4).  The  shock  layer  is  com- 

7  8  9 

puted  first  by  a  technique  which  I  developed  *  *  .  No  claim 
is  made  here  for  it  as  the  only  possible  technique  or  the  best 
one.  However,  a  certain  effort  has  been  made  to  compute 
the  boundary  points  with  a  physical  insight: 

(1)  The  pressure  is  determined  by  an  equation  where 
the  relevant  parameters  are  the  pressure  and  the  normal 
component  of  the  velocity  because  the  pressure  waves  sent 
to  the  wall  in  the  normal  directions  bounce  back  from  it  due 
to  the  condition  of  impermeability  (V  =  o).  The  technique 
used  is  a  modified  method  of  characteristics  in  the  plane 
defined  by  the  normal  to  the  wall  and  the  time  axis. 

(2)  The  entropy  is  kept  constant  along  the  particle 
path  on  the  wall;  from  pressure  and  entropy,  the  density 
follows . 

(3)  The  tangential  velocity  is  determined  by  a 
second-order  accurate  finite -different  formulation  of  a 
momentum  equation. 

A  7  x  12  mesh  is  used  to  cover  the  computational 

region.  The  final  results  are  reported  at  pages  48  and  49 

a 

of  the  report  mentioned  above  and  can  be  considered 
accurate  at  least  to  within  one  tenth  of  1%.  (Much  better 

9 

results  arc  obtained  with  a  10  x  14  mesh  --  See  pages 
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50  and  51).  Fig.  2  shows  the  pressure  distributions  between 
body  and  shock  along  five  different  radii  (9  =  0  is  the  sym¬ 
metry  line  and  0  increases  from  the  symmetry  line  to  the 
shoulder  of  the  body).  Here  again  the  non- symmetrical 
behavior  of  p  with  respect  to  the  rigid  wall  is  evident  (and 
is  qualitatively  similar  to  the  incompressible  case  mentioned 
above) .  It  may  be  expected  that  in  repeating  the  computation 
after  using  symmetry  rules  at  the  body  points  (and,  of  course, 
leaving  the  rest  unchanged)  the  results  should  not  be  too  much 
affected  near  the  symmetry  line  but  should  show  a  definite 
worsening  with  increasing  6.  This  is  what  actually  happens. 
Fig.  3  shows  the  same  lines  as  Fig.  2  as  dashed  lines,  and 
the  results  of  the  second  computation  as  solid  lines.  The 
slope  of  the  p- curves  at  the  wail  is  larger  than  the  "exact" 
one,  as  a  reaction  of  the  numerical  scheme  to  an  attempt 
to  enforce  a  vanishing  slope.  Wiggles  are  generated,  which 
propagate  through  the  shock  layer,  modify  the  shock  shape 
and  eventually  get  trapped  between  shock  and  body.  Fig.  3 

has  been  drawn  at  the  560th  computational  step,  which  is 

o 

the  last  step  computed  .  After  it,  the  results  are  prac¬ 
tically  frozen  in  both  computations.  In  cases  where  the  wall 
geometry  is  more  complicated,  the  errors  generated  by  the 
reflection  technique  may  eventually  make  the  results  un¬ 
stable. 
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Another  technique  for  prescribing  values  on  the  extra 
row  of  points  consists  of  extrapolating  from  inner  points. 
Regardless  of  the  order  of  curve-fitting  used  to  extrapolate, 
it  is  easy  to  see  that  such  a  procedure  is  arbitrary,  with  no 
connection  with  the  physical  behavior  of  the  flow.  Consider 
three  points,  A,  B,  and  C.  The  first  is  an  interior  point, 
the  second  is  the  point  on  the  wall,  and  the  third  belongs 
to  the  auxiliary  row,  behind  the  wall.  We  know  the  values 
of  a  physical  parameter,  f  at  A  and  B  and  its  value  at  C  is 
obtained  by  extrapolation.  In  going  from  t  to  the  next  com¬ 
putational  line,  t  +  A  t,  some  derivatives  are  needed,  for 
example  df/dn,  daf/dna  where  n  is  a  coordinate  in  the  nor¬ 
mal  direction.  The  increment  in  f  at  B  depends  on  those 
derivatives,  which  in  turn  depend  on  the  geometrical  nature 
of  the  fitting  curve,  not  on  physical  properties. 

For  the  sake  of  brevity,  I  do  not  present  here  any 
result  obtained  by  using  extrapolation  techniques  instead 
of  the  technique  mentioned  above  .  They  are  particularly 
bad  when  strong  curvature  effects  dominate  the  behavior  of 
the  flow  in  the  vicinity  of  the  wall.  A  more  complete  dis¬ 
cussion  of  such  effects  will  be  published  elsewhere. 

Other  attempts  to  deal  with  boundary  points  on  a  rigid 
wall,  where  the  physical  nature  of  the  problem  is  not  satisfied, 
can  be  found  in  the  literature.  Again  in  the  blunt  body  problem, 
Lapidus^  introduces  a  diffusion  along  the  boundary  by  averaging 
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certain  values  and  arbitrarily  corrects  the  momentun  vector 
to  make  it  parallel  to  the  wall.  The  results,  as  the  paper 
itself  shows,  are  not  very  good. 

3.  PERMEABLE  SUBSONIC  BOUNDARIES 

The  second  problem  to  be  carefully  considered  refers 
to  the  handling  of  those  boundaries  which  are  not  rigid  walls 
but  arbitrary,  permeable  limits  of  the  computational  region. 
For  the  sake  of  an  argument,  let  us  confine  ourselves  to  the 
case  of  an  inviscid,  compressible  flow. 

The  computational  region  must  be  finite.  Since  most 
of  the  practical  problems  deal  with  flows  in  an  infinite  domain, 
one  is  attempted  to  consider  only  a  part  of  it  by  drawing  ar¬ 
bitrary  boundaries  such  as  the  ones  shown  in  Fig.  4.  The 
flow  may  enter  the  computational  region  or  exit  from  it  through 
8 uch  boundaries.  In  Fig.  4,  CD,  GH  are  "entry"  boundaries, 
whereas  AB,  EF,  and  HL  are  "exit"  boundaries. 

If  the  flow  is  supersonic  at  all  points  of  one  of  these 
lines,  there  are  no  difficulties.  If  the  line  is  an  entry  bound¬ 
ary,  all  values  can  be  prescribed  along  it.  Making  them  con¬ 
sistent  with  the  flow  upstream  of  the  boundary  is  a  task  inde¬ 
pendent  of  the  numerical  computation.  Once  such  values  are 
prescribed,  there  is  no  feedback  into  the  region  upstream 
from  the  computational  region.  If  the  line  is  an  exit  boundary, 
there  is  no  feedback  into  the  computational  region  from  the 
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region  downstream  and  there  is  no  feedback  either  from 
the  exit  boundary  onto  the  preceding  row  of  points,  so  that  the 
values  at  the  exit  boundary  points  can  be  determined,  at  each 
step,  by  a  simple  linear  extrapolation  from  the  inside  of  the 
computational  region. 

If  the  flow  is  subsonic  on  an  entry  line,  any  change  in 
the  computational  region  will  send  waves  upstream  through  the 
line.  These  waves  modify  the  upstream  flow  and  the  values 
on  the  entry  line.  At  the  next  computational  step,  one  should 
assume  updated  values  on  that  line,  but  this  implies  bringing 
in  the  influence  of  the  upstream  flow.  Thus,  at  least  a  part 
of  the  outer  flow  should  be  treated  in  the  same  way  as  the 
inner  flow,  which  is  inconsistent  with  the  existence  of  a 
partition.  Similar  considerations  can  be  made  for  a  sub¬ 
sonic  exit  line. 

Any  permeable  limit  of  a  computational  region  with 
subsonic  points  and  constant  values  assumed  on  it  makes  a 
time-dependent  problem  ill-posed.  Disturbances  are  created 
which  propagate  inwards.  Here  again,  there  are  cases  where 
the  results  do  not  look  so  bad,  because  the  arbitrary  boundary 
is  far  from  the  region  of  interest  and  the  computation  is 
halted  before  the  perturbation  waves  have  a  chance  of  piling 
up.  But  this  is  not  a  valid  argument. 

The  only  possible  way  out  is  to  extend  the  computational 
region  to  where  the  physical  values  are  well  defined  and  un¬ 
affected  by  traveling  waves,  either  going  inwards  or  outwards. 
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This  amounts  to  extending  the  region  to  infinity.  In  certain 
cases  it  is  possible  to  do  it.  For  example,  in  the  case  of 
a  nozzle  with  a  subsonic  entrance,  one  may  imagine  the 
nozzle  extended  to  an  infinite  reservoir,  where  the  velocity 
vanishes  and  the  other  physical  parameters  take  on  their 
stagnation  values.  In  the  case  of  a  flow  about  an  obstacle, 
subsonic  and  uniform  at  infinity,  the  computation  must  be 
performed  on  the  whole  space.  We  can  do  this  and  yet  keep 
the  number  of  mesh  points  within  reasonable  limits  by  mapping 
the  infinite  physical  domain  into  a  finite  region,  a  rectangle, 
say.  Suppose  that  AB  (Fig.  5)  maps  the  infinity  of  the  physical 
plane,  where  all  physical  parameters  are  given,  constant 
values.  The  CD  line  maps  some  line  lying  in  the  physical 
plane  at  a  large,  but  finite  distance  from  the  region  of  in¬ 
terest.  The  computation  is  performed  on  the  rectangular 
mesh,  obviously  using  the  equations  of  motion  after  trans¬ 
formation  in  the  auxiliary  plane.  To  compute  the  values  on 
the  CD  line  we  use  information  from  the  AB  line  and  from 
the  EF  line.  The  information  from  the  AB  line  brings  in  exact 
values.  The  geometrical  stretching  between  CD  and  AB  (and, 
consequently,  the  damping  of  outgoing  waves  and  the  absence 
of  incoming  waves)  is  included  in  the  coefficients  which  affect 
the  transformed  equations.  One  may  be  skeptical  about  the 
amount  of  error  introduced  by  a  stretching  where  a  local 
value  of  a  coefficient  should  accurately  represent  the  global 
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phenomenon  actually  taking  place  throughout  an  infinite  region 
(such  as  the  physical  region  mapped  onto  the  strip,  ABCD). 
However,  the  flow  in  the  strip  is  nearly  uniform.  All  de¬ 
rivatives  are  small  and  the  local  evaluation  of  the  coefficients 
is  sufficient  to  insure  accuracy.  The  important  point  is  that 
the  calculation  in  the  rectangle  is  physically  well-posed. 

1  am  indebted  to  Professor  M.  Van  Dyke  for  pointing 
out  to  me  after  the  oral  presentation  of  this  paper  that  a 
similar  conclusion  was  reached  by  Wang  and  Longwell**. 
There  the  authors  compare  the  results  of  two  calculations  of 
incompressible  viscous  flow,  one  with  an  arbitrary  trunca¬ 
tion  of  the  computational  region  and  the  other  with  a  stretch¬ 
ing  of  coordinates  as  described  above.  The  results  of  the 
latter  computation  are  undoubtedly  closer  to  physical  reality 
than  those  of  the  former  one.  Because  of  the  nature  of  the 
problem,  both  sets  of  results  are  smooth.  In  the  first  case 
such  smoothness  provides  a  good  example  of  how  dangerous 
it  may  be  to  infer  accuracy  from  smoothness. 

For  compressible  inviscid  flows  wiggles  build  up  as 
a  consequence  of  the  presence  of  subsonic  boundaries.  For 
example,  Fig.  6  shows  the  pressure  distribution  along  the 
wall  of  a  two-dimensional  Laval  nozzle  as  assumed  initially 
(dotted  line)  and  after  300  computational  steps.  The  tech- 

9 

nique  used  is  the  same  as  the  one  outlined  for  the  blunt  body 
problem.  The  wiggles  in  the  results  are  a  direct  conse- 
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quence  of  an  arbitrary  limitation  of  the  computed  region.  Waves 

traveling  upstream  are  reflected  at  the  boundary  and  trapped 

within  the  computational  mesh.  Results  for  an  axisymmetric 

nozzle  computed  with  the  same  technique  but  with  a  stretching 

of  coordinates  which  let  the  nozzle  emerge  from  an  infinite  reser- 

12 

voir  are  reported  by  Migdal  et  al.  and  appear  to  be  extremely 
good. 

Another  test  of  the  infinite- stretching  idea  is  currently 

in  progress.  The  two-dimensional  compressible  inviscid  flow 

about  a  circle  in  a  stream,  uniform  at  infinity,  is  being  studied. 

For  very  low  values  of  the  free  stream  Mach  number  (M  =  0.1, 

say),  the  flow  is  practically  incompressible.  The  problem  is 

analyzed  as  the  time -dependent  evolution  of  a  flow,  starting 

impulsively  from  rest.  At  the  rigid  boundary,  the  conditions 

are  treated  as  outlined  in  the  preceding  section.  A  second- 

order  accurate  technique  is  used  for  the  interior  points.  A 

number  of  stretching  functions  has  been  used  in  an  attempt 

1 3 

to  optimize  the  clustering  of  mesh  points  near  the  body. 

The  dots  in  Fig.  1  of  this  paper  show  the  computed  values 
of  pressure  and  velocity  at  9  =  tt/2,  after  500  computational 
steps.  The  agreement  with  the  theoretical  values  for  an  in¬ 
compressible  flow  is  extremely  good,  which  proves  that  both 
boundary  conditions,  at  the  rigid  wall  and  at  infinity,  are 
properly  treated. 
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4.  IMBEDDED  SHOCKS 


A  third  problem  which  should  be  mentioned  when 
talking  about  the  treatment  of  boundary  conditions  is  the 
problem  of  discontinuity  surfaces.  I  will  confine  myself 
to  a  brief  analysis  of  shock  waves.  A  shock  can  be  con¬ 
sidered  as  a  boundary  between  two  regions  where  the  physical 
parameters  are  continuous  and  differentiable.  In  inviscid 
flows,  shocks  are  discontinuities  and  should  be  handled  as 
such  since  the  jump  conditions  across  a  shock  are  not 
compatible  with  Euler's  differential  equations.  One  may 
anticipate,  thus,  that  any  attempt  to  compute  a  flow  con¬ 
taining  a  shock  by  a  numerical  technique  which  closely 
approximates  the  inviscid  differential  equations  is  bound 

to  fail  (Richtmyer  and  Morton1,  pages  329,  Fig.  12.4; 

2 

Moretti  ).  It  is  well-known,  however,  that  sharp  trans¬ 
itions  leading  from  the  state  ahead  of  a  shock  to  the  state 
behind  it  can  be  obtained  as  solutions  of  the  Navier-Stokes 
equations  and  that  the  thickness  of  the  transitional  region 
is  related  to  the  Reynolds  number  of  the  flow.  A  similar 
result  is  obtainable  by  using  difference  equations  obtained 
from  the  partial  differential  equations  for  inviscid  flow, 
with  the  additions  of  an  artificial  viscosity  (Richtmyer 
and  Morton*,  page  327,  Fig.  12,3).  However,  the  Rey¬ 
nolds  number  is  so  high  in  most  practical  cases  that,  to  ob¬ 
tain  a  good  value  of  the  shock  thickness,  one  should  use  mesh 
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sizes  many  orders  of  magnitude  smaller  than  the  mesh  sizes 
permitted  by  the  capacity  and  speed  of  our  present  computers. 
A  coarse  mesh  and  a  high  viscosity  (either  natural  or  artificial) 
spreads  the  shock  thickness  unrealistically.  In  addition,  the 
high  viscosity  affects  the  whole  flow  field  and  makes  the  sol¬ 
ution  depart  substantially  from  the  inviscid  one  (or  the  one 
corresponding  to  the  actual,  high  Reynolds  number) .  In  this 
connection,  a  note  of  warning  should  be  given  to  the  reader 
of  page  328  in  Richtmyer's  book.  There  the  problem  is  that 
of  a  jump  between  two  constant  states,  where  all  derivatives 
vanish.  Obviously,  the  presence  of  viscosity  does  not  affect 
a  uniform  flow  and  does  not  affect  a  numerical  computation 
of  a  uniform  flow  either.  However,  the  numerical  solution 
of  any  other  case  treated  as  an  inviscid  flow  with  the  addi¬ 
tion  of  artificial  viscosity  is,  in  all  practical  calculation, 
extremely  poor  from  the  point  of  view  of  accuracy. 

The  Lax-Wendroff  scheme  (Richtmyer  and  Morton*, 
page  330)  is  undoubtedly  a  closer  approximation  to  the 
differential  equations  of  inviscid  flow  than  a  first-order 
scheme  with  added  artificial  visco  lity.  In  addition,  it  is 
consistent  with  the  Rankine-Hugoniot  conditions  across  a 
shock.  Therefore,  it  is  considered  a  powerful  tool  for  the 
numerical  analysis  of  shocked  flows,  and  some  interesting 
results  have  indeed  been  obtained.  However,  aside  from  the 
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unrealistic  spreading  of  the  shock  thickness  over  several 
mesh  points,  it  shows  another  disturbing  feature,  that  is, 
at  least  on  one  side  of  the  shock  there  are  oscillations. 

Such  a  situation  is  clearly  shown  by  Fig.  12.  6  of  Richt- 
myer's  book*,  page  333.  A  reader,  familiar  with  the 
theory  of  Fourier  series,  is  reminded  of  a  similar  pheno¬ 
menon  which  appears  any  time  a  Fourier  expansion  is 
attempted  on  a  discontinuous  function.  See,  for  example, 
Fig.  7  where  the  Fourier  expansions  of  a  saw-tooth  function, 
truncated  at  the  tenth  term  and  at  the  50th  term,  are  plotted 
in  the  vicinity  of  the  discontinuity.  The  fact  that  such  a 
Fourier  expansion  is  not  a  suitable  way  of  handling  a  dis- 

14 

continuity  is  proved  by  the  Gibbs  phenomenon  . 

Now,  when  a  Fourier  expansion  of  a  discontinuous 
function  is  needed,  one  proceeds  first  to  eliminate  the  dis¬ 
continuities  by  subtracting  from  the  given  function  the  pro¬ 
per  number  of  saw-tooth  functions  with  suitable  values  of 
the  jumps,  and  then  to  expand  the  continuous  function  so  ob¬ 
tained.  The  convergence  is  quick  and  the  Gibbs  phenomenon 
does  not  appear. 

The  same  principle  should  inspire  the  numerical 
treatment  of  fluid  flows  containing  shocks.  The  regions 
between  shocks  are  smooth  and  can  be  handled  with  coarse 
meshes.  The  shocks  (which,  in  a  time-dependent  case, 
must  be  left  free  to  move)  can  be  described  with  ease  by 
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the  Rankine-Hugoniot  conditions  plus  one  equation  written 

from  the  inside  of  each  region  up  to  the  shock.  Such  an 

equation  must  have  the  nature  of  the  compatibility  equation 

on  a  characteristic.  In  multi -dimensional  problems,  the 

method  of  characteristics  has  to  be  modified  suitably.  That 

9 

this  can  be  done  is  proved  by  the  excellent  results  obtained 
if  all  the  shock  and  body  points  are  computed  using  such  a 
technique . 

15 

Last  year,  1  reported  a  preliminary  approach  to  a 
more  ambitious  problem.  A  steady  three-dimensional  sym¬ 
metric  flow  was  studied,  which  is  strongly  compressed  in  a 
region  and  strongly  expanded  in  another  one.  An  imbedded 
shock  builds  up  in  a  part  of  tire  flow  and  dies  out  somewhere. 

in  a  simpler  case,  the  possibility  has  been  proved 
not  only  of  fitting  a  moving  shock  so  that  the  solution  is 
extremely  accurate  despite  the  use  of  a  coarse  mesh,  but 

also  of  predicting  where  and  when  the  shock  starts  building 

2 

up  without  resorting  to  coalescence  oi  characteristics  .  in 
Fig.  3,  some  frames  from  a  computer-made  movie  are  shown, 
from  which  the  formation  of  a  shock  appears  together  with  its 
evolution  and  final  freezing  at  the  place  where  the  steady  shock 
should  be,  according  to  the  geometry  of  the  nozzle  and  the 
exit  pressure. 
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5.  CONCLUSIONS 


Let  us  summarize  now  the  basic  points  of  the  dis¬ 
cussion. 

(1)  There  are  two  classes  of  points  in  a  numerical 
computation,  interior  points  and  boundary  points. 

(2)  No  conclusion  can  be  drawn  on  the  effectiveness 
of  a  technique  for  interior  points  unless  one  is  sure  that  no 
disturbances  are  fed  in  from  the  boundaries. 

(3)  One  can  discuss  the  degree  of  accuracy  of  a 
particular  technique  for  interior  points;  however,  within 

the  current  state  of  the  art,  the  question,  as  fa.r  as  boundary 
points  are  concerned,  is  different,  viz.;  Is  the  technique 
used  right  or  wrong? 

(4)  Cures  for  "non-linear  instability"  such  as  artificial 
viscosity  decrease  accuracy  and  conceal  the  nature  of  the 
trouble. 

(5)  Rigid  wall  points  can  be  properly  treated.  Some 
examples  are  shown  in  this  paper  and  warrant  more  intensive 
study. 

(6)  Preliminary  examples  show  that,  for  certain  flows 
with  simple  conditions  at  infinity,  the  use  of  a  stretching  tech¬ 
nique  provides  satisfactory  results. 

(7)  Imbedded  shock  waves  should  be  treated  as  dis¬ 
continuities,  if  one  aims  to  achieve  accuracy  without  having 
to  resort  to  a  prohibitively  fine  computational  mesh. 

(8)  Beautiful  results  can  be  achieved  on  the  present 
generation  of  computers,  using  a  small  number  of  mesh 
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points  and  consequently  spending  a  short  time  to  compute,  if 
the  boundary  value  problem  is  properly  handled.  In  fact,  I 
believe  that  finite -difference  techniques  for  initial-and- 
boundary-value  problems  are  a  powerful,  accurate,  and 
inexpensive  computation  tool. 


ACKNOWLEDGMENT 

This  research  was  conducted  in  part  under  Contract 
Nonr  8  39(34),  and  under  Contract  Nonr  839(38)  for  PROJECT 
STRATEGIC  TECHNOLOGY,  and  was  made  possible  by  the 
support  of  the  Advanced  Research  Projects  Agency  under 
Order  Nt  529  through  the  Office  of  Naval  Research. 


22 


REFERENCES 


1.  Richtmyer,  R.  D.  and  Morton,  K.  W.,  Difference 
methods  for  initial  value  problems,  Second  Edition, 
Interscience  Publ. ,  New  York,  1967. 

2.  Moretti,  G.  ,  Starting  a  Laval  nozzle  --  A  study  of 
numerical  techniques,  to  be  published  as  a  Polytechnic 
Institute  of  Brooklyn  Report. 

3.  Bohachevsky,  I.  O.  and  Rubin,  E.  L.  ,  A  direct  method 
for  Computation  of  Non-equilibrium  flows  with  detached 
shockwaves,  AIAA  J.  ,  4,  600,  1966. 

4.  Bohachevsky,  I.  O.  and  Mates,  R.  E.,  A  direct  method 
for  calculation  of  the  flow  about  an  axisymmetric  blunt 
body  at  angle  of  attack,  AIAA  J. ,  4,  776,  1966. 

5.  Burstein,  S.  Z.  ,  Numerical  Methods  in  multi-dimensional 
shocked  flows,  AIAA  J.,  2,  2111,  1964. 

6.  Singleton,  R.  E.  ,  Lax-Wendroff  difference  scheme  applied 
to  the  transonic  airfoil  problem,  AGARD  Specialists' 

Meeting  on  Transonic  Aerodynamics,  AGARD  Con¬ 
ference  Proceedings  No.  35  (Reprint),  18-20  September 
1968. 

7.  Moretti,  G.  ,  and  Abbett,  M.  ,  A  time -dependent  com¬ 
putational  method  for  blunt  body  flows,  AIAA  J.  4,  2136,  1966. 

8.  Moretti,  G.  and  Bleich,  G. ,  Three-dimensional  flow  around 
blunt  bodies,  AIAA  J. ,  5,  1557,  1967. 

9.  Moretti,  G.  ,  Inviscid  blunt  body  shock  layers  --  Two- 
dimensional  symmetric  and  axisymmetric  flows,  PIBAL 
Report  No.  68-15,  1968. 


23 


10.  Lapidus,  A. ,  A  detached  shock  calculation  by  second 
order  finite  differences,  J.  Comp.  Phys.  2,  154,  1967. 

11.  Wang,  Y.  L.  and  Longwell,  P.  A.,  Laminar  flow  in  the 
inlet  section  of  parallel  plates,  A.  I.  Ch.  E.  J.  10,  323, 
1964. 

12.  Migdal,  D.  ,  Klein,  K.  and  Moretti,  G.  ,  Time  dependent 
calculations  for  transonic  nozzle  flow  (to  be  published  as 
a  technical  note  in  the  AIAA  J. )  „ 

13.  McKenzie,  D.  and  Moretti,  G. ,  Time -dependent  calcula¬ 
tion  of  the  compressible  flow  about  airfoils.  AGARD 
Specialists'  Meeting  on  Transonic  Aerodynamics,  AGARD 
Conference  Proceedings  No.  35  (Preprint),  18-20 
September  1968. 

14.  Moretti,  G. ,  Functions  of  a  complex  variable ,  Prentice- 
Hall,  New  York,  1964,  page  334. 

15.  Moretti,  G.  and  Bastianon,  R. ,  Three-Dimensional 
effects  in  intakes  and  nozzles,  Paper  No.  67-224, 

AIAA  5th  Aerospace  Science  Meeting,  New  York,  1967. 


24 


29.2° 


SHOCK 


FIG.  3 


2 


2.0 


31 


Unclassified _ 

^ — ^ ^ ^ _ _ _ 

DOCUMENT  CONTROL  DATA  .  R  &  0 


(Security  elmaalllcmtlon  ol  title,  txuly  ol  nhntnict  imd  tinh'Mt' irf  mum  tut  Inn  must  he  antetrtl  when  the  overall  report  1 1  rltt**l Itfd) 


1  ORIGINATING  ACTIVITY  fCorpormlt  ntilhtyr) 

Polytechnic  Institute  of  Brooklyn 

Dept,  of  Aerospace  Engineering  and  Applied  Mechanics 
Route  110,  Farmingdale,  New  York  11735 

i0,  REPORT  SECURITY  CLASSIFICATION 

Unclassified 

2h.  GROUP 

J  REPORT  TITLE 

THE  IMPORTANCE  OF  BOUNDARY  CONDITIONS  IN  THE  NUMERICAL  TREATMENT 
OF  HYPERBOLIC  EQUATIONS 

A  OE5CRIPTIVE  NOTES  (T\po  ol  repot t  nnd  Inclusive  ttntee) 

Research  Report 

3  A u  T ho R  1*1  (  b'lrat  nemo,  middle  Initial,  laat  t»amv) 

Gino  Moretti 

6  REPOR  T  O A  T  C 

November  1968 

7  n.  total  no  oi  pages  7fi.  no  o  r  n  r  r ' 

34  15 

B«.  CONTRACT  OR  GRANT  NO 

Nonr  8  39(34)  and  Nonr  839(38) 

6.  PROJEC  T  NO 

NR  061-135 

c.  ARPA  Order  No.  529 

d. 

*J#».  ORIGINATOR'S  ME  PONT  HUM  HI  Rl5> 

PIBAL  Report  No.  68-34 

•»h.  OTHER  R». port  NOISI  (Any  other  number*  thot  mm-  !»••  U'stftnuJ 
thla  report) 

10  DISTRIBUTION  STATEMENT 

Distribution  of  this  document  is  unlimited. 

II  SUPPLEMENT  ARY  NOTES 

M  SPONSORING  MILITARY  ACTIVITY 

Office  of  Nava)  Research 

Department  of  the  Navy 

Washington,  D.C.  20360 

1  13  ABSTRACT  1 

*Many  of  the  existing  computations  of  initial- and -boundary  value  problems  in  fluid 
mechanics  suffer  from  unrealistic  treatment  of  boundary  points.  Three  categories 
of  boundaries  are  briefly  discussed:  rigid  walls,  arbitrary  boundaries  of  a  computa¬ 
tional  region  in  a  subsonic  flow,  and  shock  waves.  An  attempt  is  made  to  show  in  what 
sense  the  numerical  treatment  of  such  boundaries  may  be  physically  wrong  and  what 
can  be  done  instead.  Examples  from  the  blunt  body  problem,  the  transonice  flow  in  a 
nozzle,  the  incompressible  inviscid  flow  past  a  circle,  and  the  quasi-one-dimensional 
flow  in  a  Laval  nozzle,  are  shown. 


DD  ,”“.1473 


Unclassified 

Sri  <lfl1\  t'l.i'.siln  ..  1  mu 


Unclassified 

S*Turny  Cl.ts-.ifi.. ill. >11 


